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Abstract. A stochastic model of the chemical enrichment of metal-poor systems by core-collapse (Type II) supemovae is 
presented, allowing for large-scale mixing of the enriched material by turbulent motions and cloud collisions in the interstellar 
medium. Infall of pristine material is taken into account by following the evolution of the gas density in the medium. Analytical 
expressions were derived for the number of stars enriched by a given number of supernovae, as well as for the amount of mass 
with which the ejected material from a supernova is mixed before being locked up in a subsequently formed star. It is shown 
that for reasonable values of the gas density (~ 0.1 cm -3 ) and of the supernova rate (~ 0.25 kpc -3 Myr -1 ) of the Galactic 
halo, the resulting metallicity distributions of the extreme Population II stars show a distinct cut-off at [Fe/H] ~ —4. In fact, 
by assuming no low-mass Population III stars were able to form out of the primordial interstellar medium, the derived fraction 
of stars below [Fe/H] = —4 is in agreement with observations. Moreover, the probability is high that even the most metal-poor 
stars observed to date have been enriched by several contributing supernovae. This partly explains the relatively small star- 
to-star scatter in many chemical-abundance ratios for stars down to [Fe/H] = —4, as recently found in several observational 
studies. Contribution from the thermonuclear (Type la) supernovae is found to be negligible over almost the entire extremely 
metal-poor regime. Although the fraction of contaminated stars may increase rapidly towards [Fe/H] = —2.5, the fraction of 
stars with iron primarily from Type la supernovae remains small. The stars that are heavily polluted by Type la supernovae are 
pushed towards higher metallicities, creating a hole/bump in the metallicity distribution. Such features could be used to reveal 
the possible presence of subpopulations of Galactic halo stars that have been enriched by Type la supernovae. 

Key words. Stars: abundances - Stars: Population II - Stars: statistics - ISM: supernova remnants - Galaxy: halo - Galaxies: 
evolution 



£> . 1. Introduction Could it be that averaging effects were already in play at 

a metallicity of [Fe/H] ~ —4? This question leads us to an- 
It is an interesting fact that in the extremely metal-poor regime other: What are the number and the metallicities of stars en- 
; ([Fe/H] < -2.5) of the Galactic halo, the observed star-to- riched by any given number of SNe? Several chemical evolu- 
star scatter in abundance ratios of elements heavier than he- tion models have been developed to study the effect of small 
lium tend to be smaller than previously anticipated (Carretta number statistics on the build-up of chemical elements in the 
et al. 2002; Cayrel et al. 2004), while in the metal-rich regime early Galaxy (e.g., Tsujimoto et al. 1999; Argast et al. 2000; 
([Fe/H] > -1) of the Galactic disk, the scatter seems larger Oey 2000; Travaglio et al. 2001). However, these models do 
than one would expect from simple stochastic enrichment con- not explicitly trace or discuss the metallicity distributions of in- 
siderations (e.g., Edmunds 1975). As noted by, e.g., Cayrel et dividual populations of stars enriched by a specific number of 
al. (2004), the dispersion is surprisingly small in many abun- SNe. In the simplest case one might assume that the metallic- 
dance ratios, excluding the neutron-capture elements, even for ity is proportional to the number of enriching SNe (e.g., Fields 
stars with a metallicity as low as [Fe/H] = -4. Supposedly e t al. 2002); but as we shall see, allowing a variable dilution 
(see, e.g., Woosley & Weaver 1995), the stellar yields of differ- (increasing with time) of the SN material and a variable stellar 
ent heavy nuclei depend strongly, among other things, on the yield introduces a dispersion in the number of SNe contributing 
progenitor mass of the supernova (SN). It therefore seems that heavy elements to stars, even for stars with identical metallic- 
a small scatter can only be obtained if these stars were formed ity. Moreover, due to the continuous mixing of the interstellar 
out of gas enriched by several SNe (cf. Karlsson & Gustafsson medium (ISM) and the possible infall of primordial gas, stars 
2005, hereafter Paper II). of a given metallicity will, on average, be enriched by more 

SNe than expected from such simple considerations. 
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The interesting approach taken by Oey (2000) examined 
stochastic enrichment of the ISM by randomly distributed star- 
forming regions, where each region was assigned a specific 
metallicity depending on its predetermined size. A similar ap- 
proach was used here. However, we focus instead on the SN 
remnants themselves, which allows us to study the effects of in- 
dividual stars on the chemical enrichment of the ISM. A more 
detailed comparison between the present approach and that of 
Oey (2000) will be given in Sect. 151)1 

In Sect. |3 the stochastic theory is outlined, including a 
derivation of the expression for the number of stars enriched 
by a certain number of SNe, as well as an expression for 
the amount of dilution of the newly synthesized material. A 
stochastic enrichment model for metal-poor systems is then 
presented in Sect. [5] making use of these results. In Sect. |4] 
three specific models are presented, applicable to the Galactic 
halo, and in Sect. [5] the results are discussed. A summary is 
given in Sect. [6] 



2. The chemical evolution of the ISM 

The chemical enrichment of the ISM proceeds as the material 
ejected in SN explosions, planetary nebulae, and stellar winds 
disperses and gradually mixes with the ambient medium. Due 
to the non-instantaneous mixing, local enrichment events, such 
as SNe, cause the interstellar abundance A/H of an element A 
relative to hydrogen to vary over space and time. This variation 
can be described by a function /a(£, x), which is different for 
each element. In order to link the evolution of chemical inho- 
mogeneities in the ISM with the abundances observed in low- 
mass stars, we require the following knowledge: what and how 
much of an element is ejected from each SN; when, where, 
and at what rate SNe explode; how the material is dispersed; 
and finally, when and where the local abundance is probed by 
subsequent star formation. Mapping of x) onto the popu- 
lation of low-mass stars is controlled and weighted by the star 
formation rate (SFR) ip = ip(t, x), which is thought to favour 
cool, dense, and gravitationally unstable regions of the ISM. 
This last step contains an integration over space and time and 
results in abundance distributions describing the total probabil- 
ity of finding a star with a certain set of chemical abundances. 
Such distributions can be compared with the observations. 

Modelling the evolution of chemical elements with intrin- 
sic fluctuations mapped by the low-mass stars is an intricate 
task. Even extensive numerical simulations suffer from sim- 
plified physics and inadequate spatial resolution. In this study, 
we instead discuss an analytical approach that neglects the de- 
tailed physics of star formation and treats the discrete chemical 
enrichment as a stochastic process. 

2. 1. Properties of the model ISM 

Let us first write down some definitions. The probability that a 
newly formed star has a mass m is given by the initial mass 
function (IMF), which is assumed to be a time-independent 



random variable with probability density function f m = <f>(m), 
normalized by number such that 

m u 

J (j)(m)dm = 1. (1) 

mi 

The lower and upper mass limits are denoted by mi and m u , 
respectively. 

In this study, we shall mainly discuss the early chemical 
evolution of primary elements — like oxygen, magnesium, and 
iron — produced in core-collapse SNe (i.e., mainly Type II SNe; 
however, see Sect. 15.51 for a discussion of the effect of ther- 
monuclear SNe). The total mass of such an element A, pro- 
duced in the progenitor star and in the explosion, is predomi- 
nantly determined by the mass of the star and shall be denoted 
by pa{itl). The mass fraction of A in the gas from which the 
star was formed should be taken into account formally when 
calculating the total amount of mass ca ejected in the explo- 
sion. However, at low metallicities, say for [Fe/H] < —2, the 
ejected amount of A is completely dominated by the yield pa- 
Thus, we will assume that = pa{tu). 

Star formation is a clustered phenomenon, and many stars 
are formed in groups like globular clusters, open clusters, and 
OB associations. However, there is also less clustered star for- 
mation, and most of the field stars cannot have been formed in 
strong gravitationally bound associations. A fundamental as- 
sumption in the present model is that star formation is unclus- 
tered and that the stars are randomly distributed in space. This 
means that the SFR can be described by a function t/> = ip(t) 
that depends only on time. Under this assumption, chemical en- 
richment events, such as SNe, are spatially uncorrelated, as is 
mapping of the chemical inhomogeneities. Within the present 
theory this assumption may, however, be replaced by the as- 
sumption of correlated star formation in stochastically indepen- 
dent star-forming regions (Karlsson 2005, hereafter Paper III). 

Like the SFR, the density p(t, x) in the model ISM is de- 
scribed by the spatial average p(t). The density changes over 
time due to three possible mechanisms: infall of gas onto the 
system, outflow of gas, and star formation. The rate of change 
is governed by 

^E. - - - — n\ 

dt ~ Uacc Uc) dt ' K ) 

where u acc is the accretion rate, u C j the ejection rate, and ds /dt 
the rate of mass locked up in stars and stellar remnants, which 
is given by the amount of mass in newly formed stars minus the 
amount of mass returned to the ISM per unit time by the dying 
stars. The lifetime r of a star depends mainly on its mass and 
will therefore be denoted by the function r = g T (m). Using 
this notation, the stars of mass m returning their outer layers 
to the ISM at time t were formed at a rate ijj(t — g T (m)) per 
unit volume. Moreover, the mass of the least-massive dying star 
depends on time such that mt(t) = Hence, 

ds f 

— = (m)i/)(t) — / (m — m r )4>(m)i^(t — g T (m))dm, (3) 
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Fig. 1. Illustration of the gradual increase in the volume con- 
taining material from SNe. The time-like parameter /i, defined 
in Eq. (|5}, is a measure of the average number of SNe that have 
enriched a volume element. Here, /i a f 2 . The SNe are ran- 
domly distributed over space at a constant rate us~n- The shaded 
bar indicates the number of enriching SNe. Time is expressed 
in arbitrary units. 



where (to) = f m u m<f)(m)dm is the average mass of stars in 
a stellar generation. Note that mass m r of the stellar remnant 
is also a function of m. Given the infall- and outflow rates, to- 
gether with an initial value po, p(t) is fixed by Eq. Q. The 
gas density will be used in calculating the amount of dilution 
of the newly synthesized elements (ejected in the SN explo- 
sions) before they are locked up in subsequently formed stars 
(Sect.E3J. 

The way infall is introduced in Eq. (|2ji implies that the mix- 
ing of infalling gas with the ISM is instantaneous. The infall is 
assumed to be a homogeneous, diffuse flow of gas onto the sys- 
tem and not in the form of, e.g., discrete clouds. We will also 
suppose that the infalling gas is pristine, consisting of pure hy- 
drogen and helium, which implies that any chemical inhomo- 
geneities produced by our model originate and evolve within 
the system. Thus, evolution beyond the system is not consid- 
ered. 

2.2. The mixing volume picture 

The total volume enriched by debris from a single SN shall 
be defined as the mixing volume V m i x of that SN. The mixing 
volume will grow steadily with time, initially due to expan- 
sion of the SN remnant and later on due to large-scale mixing 
processes, such as turbulent diffusion. In this picture a mix- 
ing volume is generated by each SN event. Chemical inho- 
mogeneities develop with subsequent generations of SNe that 
explode within earlier generations of mixing volumes, so that 
mixing volumes begin to overlap and the gas becomes further 
enriched. This scenario is illustrated in Fig. ^ At any instant 




& 0.05 



o.oo 

16 



16.0 



k (no. enriching SNe) 



0.0 



Fig. 2. The two-dimensional probability density function 
wism (k, t) as given by Eq. (|4}. In particular, the thick line gives 
the probability of finding a region in space enriched by 2 SNe 
as a function of the time-like variable The dots give the prob- 
ability of finding a region enriched by k SNe at fi = 2.0 (cf. 
Fig. 0. Note that u>ism(0, t) is the probability of finding an 
unenriched region. 



in time, different regions of the ISM will be enriched by dif- 
ferent numbers of SNe, and a distribution u>isM(fc, t) exists to 
describe the probability of finding k overlapping mixing vol- 
umes at time t somewhere in the ISM. Knowing this distribu- 
tion, we are able to quantify and map the abundance distribu- 
tions in stars. 

2.3. The distribution of chemical inhomogeneities 

To find a general expression for u>isM(fc, t) for a turbulent and 
inhomogeneous ISM with a spatially clustered SN distribution 
is beyond the scope of this study. However, if the SNe are ran- 
domly distributed over space with an average rate usn per unit 
time and unit volume, we take the probability of finding a re- 
gion enriched by k SNe at time t to be approximated by the 
Poisson distribution (see Fig.|2}, i.e., 



wism(M) = P(k,p,{t)) = 



p(t) k /kl 



(4) 



where /j,(t) is the integrated volume affected by SNe, expressed 
in units of the total volume of the system at time t. Dispersive 
processes, like turbulent diffusion, will gradually spread the SN 
material over larger and larger volumes. A mixing volume that 
grows with time would mimic such an effect. More precisely, 
if V m i x (t — t') is the size of the volumes at time t enriched by 
material from SNe that exploded at time t', then 



Hit) = / V mi x(t - t')u SN (t')dt'. 



(5) 



With this definition, p(t) is a measure of the average number 
of SNe that have enriched a random volume element in space 
at time t. By introducing a function of the form 6{t — t')0(t'), 
where 9 is the Heaviside step function, we can formally extend 
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the limits of integration to ±00, thus expressing /iasa convo- 
lution such that 



m(*) = / tfnu(t - t'W - t')e{t')um{t'W 

= V mix * wsn(*)j 



(6) 



where V ro i x = V m i X (t)6(t) and Usn = usu{t)0(t). This result 
will be used in Sect. 12.41 to derive an expression for the total 
number density of stars enriched by a certain number k of SNe. 



Eq. 0. The derivative of a convolution is formally given by 
^(V mix *w SN ) = V mix *u S N (= V mix *w S N) where the deriva- 
tive Of Vmiv is 



mix 

d 



dt 



(Vmi x 9) = V m i x # + V m i x 5, 



(9) 



and 5 is the Dirac delta function. Since V m i x (0) = 0, the second 
term in Eq. (|9} vanishes in the integration, and we obtain 



dfi 
~dl 



Vmix *USN = / (f - t')u SN (t')dt'. 



(10) 



2.4. Star counts 

Even though the increase in the average abundance of an el- 
ement in the ISM can be considered as a continuous evolu- 
tion, the local abundance increases due to discrete enrichment 
events, such as SNe. Therefore, it might be more natural to des- 
ignate the build-up of heavy elements in the ISM as a chemical 
enrichment rather than a chemical evolution, especially in the 
extremely metal-poor regime. We emphasize this distinction in 
our treatment by integrating over time and deriving expressions 
with the number k of enriching SNe as the explicit variable. 

Let us first calculate the total number of stars enriched by a 
given number k of SNe assuming that the probability of finding 
a region enriched by k SNe at time t is given by Eq. 0. If 
is the SFR per unit time and unit volume, the number density of 
stars formed in regions enriched by k SNe in the time interval 
[t, t+dt] is u>isM(fc, t)tp(t)dt. The total numberof still-existing 
stars per unit volume formed after a long time is then given by 
the integral 



n k = a Jr (i)wi SM (A;,t)V'(0 d ^ 



(7) 



where tq — 13.5 Gyr is the age of the oldest stars (Bennett 
et al. 2003), which we assume coincides with the age of the 
Galaxy, and a* (t) is the fraction of a stellar generation formed 
at t that still exists today. The latter is given by 



a*(£) = J (j)(m)dm, 

nil 



(8) 



where (f>(m) is the IMF and <7 1 T 1 (tg — t) denotes the mass of 
the most massive, still-surviving stars formed at time t. Note, 
however, that in the early Galaxy, when most of the stars with 
[Fe/H] < — 2.5 were formed, a* changes very little. Therefore, 
in the extremely metal-poor regime we may approximate the 
integral in Eq. (|8j with a constant such that a* = slms. where 
olms is me fraction of low-mass stars with lifetimes greater 
than the age of the Galaxy. 

A closed expression for n k cannot be obtained in the gen- 
eral case. However, by the change of variables t — > /i we 
can approximate the integral in Eq. @ by a function M and 
its derivatives. To find AT, we first need an expression for 
dfi/dt, which is given by the derivative of the expression in 



The SN rate usn(*) relates to ijj{t — tsn)> where tsn is the 
typical lifetime of a high-mass star that explodes as an SN. 
However, neglecting this small time delay between the SFR 
and the SN rate and assuming a time-independent IMF, usn is 
proportional to rp such that 



(11) 



where osn is the fraction of high-mass stars in a stellar genera- 
tion. This assumption is used in all the calculations that follow, 
and for Eq. (II 01 we get 



dfi 



osn / V mix (t - t')^{t')dt' = a SN ^, 



(12) 



where the integral denoted by 1Z = lZ(fi(t)) is a measure of 
the rate of change of /1. 

Since both V m i x and usn are, by necessity, non-negative 
functions of time, /j,(t) increases monotonically and t can there- 
fore be expressed in terms of /i. Thus, the change of variables 
t — > fJ, yields d^i/dt — asN^(M) an d V 1 = VKm), so that Eq. @ 
becomes 



n k 



1 

ki 



1>/K 



osn/olms 



d/j,, 



(13) 



where /iq = m(tg). Recall that for the early Galaxy, 
a* ~ etLMS- As soon as w\sm starts to decrease, it quickly goes 
to zero due to the exponential (see Fig.|2). This means that 
for relatively small k (i.e., at low metallicities, depending on 
the value of /ig) the upper integration limit /jg can, in practice, 
be extended to +00. Identifying the function N with 



AT( M ) = 



aSN/«LMS 

we may thus rewrite Eq. Jl 31 such that 



nh 



Urn 



1 



(14) 



(15) 



i.e., rii. is equal to the expectation value of AT. In this respect, 
wism should be regarded as a weighting function and the sub- 
script k on the bracket indicates the dependence on k. An ap- 
proximate solution to this integral can be found by expanding 
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Af(p) in a Taylor series around p = k. Including the second- 
order term, we get 

n k ~Af(k)+Af'(k) + ^(k + l)M"(k), (16) 

where, e.g., H' = dM/dp. Note that if H(p) is close to linear, 
then nfc ~ AT(k) (+const). 

There is an interesting limiting case when assuming a static 
ISM and a constant mixing volume. Let V m - lx be the typical 
volume affected by each SN such that the ejected material is 
instantaneously mixed within this volume. Using Eq. il Q , the 
expression for p(t) in Eq. 10 is then reduced to 



p{t) = asN^mix J ip(t)dt , 


A closed expression for n./. can be found such that 



(17) 



" 





d/i' = a SN V mix V>(t') d *' 



1/V n 



M 

lim I / e-M'^d/x' = 1/Fmix . (18) 



osn/olms m 



The ratio is written in this form to show the connection to N in 
Eq. J14i . This result tells us that there are as many (low-mass) 
stars enriched by 1 SN as there are stars enriched by 2 SNe, by 
3 SNe, and so on; i.e., n k is independent of k, even for a varying 
SFR! This is to be expected since usn(£) = asNip(t), meaning 
that mapping the chemical inhomogeneities by low-mass star 
formation occurs at the same rate as the inhomogeneities are 
generated. It is a limiting case since any time-dependent mixing 
volume growing with time would result in a decreasing n k with 
increasing k. This result was used without proof in Karlsson & 
Gustafsson (2001, hereafter KG). 

In their work on the evolution of SN remnants, Cioffi 
et al. (1988) find that the final radius of the remnant, when 
it merges with the ambient ISM, is proportional to p~ 18 / 49 . 
The amount of swept-up mass is thus proportional to p~ 5 / 49 . 
Let us neglect this weak dependence on density and assume 
that the material from a typical SN becomes diluted within a 
constant mixing mass. The mixing volume is then given by 
Vniix(i) = M m i x / p(t), where M m - 1K is the mixing mass and 
where p(t) is determined by Eq. (|2j. Again assuming instanta- 
neous mixing within the mixing volume, we have 



p{t) = a SN M mix / ^(t')/p{t')dt 



Hence, dp/dt cx tp/p. Using Eqs. ( H4H 16> we obtain 

l/-^xmix / \ l/-^4mix / 7 \ 
n k = ; \p)k ^ ; P(k), 



osn/olms 



osn/olms 



(19) 



(20) 



by keeping only the zeroth order term. This expression differs 
somewhat from the one in Eq. Jl 8b . Suppose that the SN rate 
(oc ijj) increases with time. The evolution of the chemical inho- 
mogeneities, characterized by the average number of overlap- 
ping mixing volumes, then speeds up. This is compensated for 
by a faster mapping of the inhomogeneities via the higher SFR 
00). However, if the density of the ISM also increases with time 
(e.g., due to infall), the size of the mixing volumes decreases 
and evolution of the inhomogeneities is slowed down accord- 
ingly. Since the SFR is unaffected, more stars are therefore able 
to map this state of chemical inhomogeneity. Thus, the number 
of stars enriched by k SNe increases with increasing k, for an 
increasing density as indicated by Eq. i2Q\ . In fact, if tp oc p, 
then n% oc ip, contrary to the result in Eq. Jl 8I >. 

In the derivation of the two last results, we assumed that the 
mixing volumes have a certain, pre-determined size and that 
the SN material is instantaneously mixed within this volume. 
However, due to the continuous agitation of the ISM the mixing 
volumes will grow with time. This means, among other things, 
that there exists a distribution of mixing masses rather than a 
single constant mixing mass. 

2.5. The distribution of mixing masses 

We now address the question of how much an element ejected 
in an SN explosion is diluted before being locked up in a 
subsequently formed star. Assuming homogeneous mixing, the 
amount of dilution is proportional to the mass of the interstel- 
lar material inside the mixing volume associated with the SN. 
For a time-dependent density p(t), the total mass M inside this 
mixing volume not depends only on the size of V,nix but also on 
when V m i x was created by the SN explosion. Let ry be the time 
that has elapsed since the creation of Vi n i x . Thus, ry denotes 
the age of V m i X and is a measure of its size. If the SN exploded 
at time t — Ty , mass M is given by the two-dimensional func- 
tion gM, where 



9M(t,T V ) = \t-r v 

0, 



_ I J VmUt' - t + T V )p{t')dt', t > T V 



(21) 



t < T V . 



The mass of enriched gas only grows by mixing of the sur- 
face of VJnix with the surrounding material. The mixing mass 
can neither change nor, in particular, decrease, due to density 
changes within V m i X . The term V m i X p is therefore not included 
in the integral in Eq. (I21> . 

Now, at time t > Ty, mixing volumes of size V m i X (Tv) 
appear at a rate usn(^ — tv), which is a measure of how 
many mixing volumes there are of this size at time t. Thus, 
the probability that such a mixing volume overlaps a region 
already enriched by k — 1 SNe at time t is proportional to 
KnixOv )wisM(k — 1, t)usj<s(t — t v )■ The probability of form- 
ing a still-existing star in this region in time interval [t, t + 
dt] and age interval [tv,tv + dry] is then proportional to 
asNa i ,(t)V m i X (Tv)wisM(k - l,t)ip(t- T V )ip(t)dtdTv, where 
we have used Eq. (II It to write wsn in terms of ijj. Integrating 
this expression over the region ADm in the (t, Ty )-plane gives 
the total probability /m 4 of having a mixing mass of size M in 
regions enriched by k SNe (cf. the treatment in KG). Thus, 
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Fig. 3. The distribution of mixing masses in regions enriched 
by k = 1, 2, 10, and 100 SNe as given by Eq. d25l for a constant 

For example, for k 



SFR (see Model A in Sect. 14.4. It , For example, for k = 10 all 
mixing masses Mi in Eqs. (I26i and J27> should be picked from 
the density function fM k=10 - 



fM h = c kj a*(t)V mix (Ty)u; IS M(fc - l,t) X 

AD M 

x - rv)^(*)d*dr F , (22) 

where the SN fraction agN has been included in the normal- 
ization constant c%. The differential integration region ADm is 
defined by function gAi(t, 7y) in Eq. (I21t . such that 



AD M = {(*,7v); Af < g M (t,T V ) < M + dM}. 



(23) 



We may simplify Eq. 1221 by assuming that the density of 
the ISM is constant. The mixing mass M is then directly pro- 
portional to the mixing volume such that 



(24) 



where ~p is a typical density. Since the mass is determined by 
the size of Vi n ; x alone, the integration over ry is redundant and 
Eq. \22\ reduces to a one-dimensional integral 



f Mk = c k M J w 1SM {k - l,t)i>{t - 7v)V(t)d*, (25) 

9m(M) 

where we used Eq. MAI to write V m i x (Ty ) in terms of M, which 
can be taken outside the integral. The factors olms an d 1/p 
have been included in the normalization constant cu- Recall that 
tq denotes the age of the Galaxy. Since mixing masses of size 
M = <7m(tv) do not exist at times t < Ty, a lower integration 
limit is introduced at t = Ty, making it, as indicated, a function 
of M. 

As shown in Fig.|5J the distribution fj^ k is broader and the 
average mixing mass is larger for higher values of k, which 
is a result of the turbulent diffusion of the ISM, mimicked by 
the growing mixing volumes. The material ejected by the very 
first SNe has continuously been diluted and spread over large 
regions. Their mixing masses are therefore large. 



3. The stochastic chemical enrichment model 

Let us now develop an inhomogeneous chemical-enrichment 
model based on the mixing-volume picture of the ISM. 
Suppose that an SN explodes inside the mixing volume of 
an earlier SN. At this point the mass fraction of element A 
in the first mixing volume is pjy{mi)/Mi, where jj^fmi) is 
the yield of A from an SN of mass mi and Mi is the to- 
tal mass within the first mixing volume. The total mass of A 
within the second mixing volume is given by the yield of A 
ejected from the second SN plus the fraction of A from the 
first SN, i.e., pa ("12) + P A(jn\)M2 j M\. The mass fraction is 
thus Pa{™>2)/M% + pa(itii)/Mi. Similarly, assuming an ini- 
tially primordial ISM, the abundance of element A in regions 
enriched by k SNe is then given by the sum 



(A/H) = 



\jrn-A PA{ m i) _ 



X H /m H 



i=l 



Mi 



= Vk(A)i 



(26) 



where mi is the progenitor mass of the ith SN, and Mj is the 
mass of the mixing volume belonging to that SN. The atomic 
masses of element A and hydrogen are denoted by rriA and mn, 
respectively. In the extremely metal-poor regime the hydrogen 
mass fraction does not change drastically from its primordial 
value, and we will assume that X\\ = 0.75 (Y p = 0.25, Huey 
et al. 2004). However, for higher metallicities a first-order cor- 
rection can be introduced such that Xn = — 3.00Z + 0.75, 
where Z is the average mass fraction of metals (based on data 
taken from Cox 2000). 

For the sake of completeness we shall also write down the 
expression for an abundance ratio of two heavy elements A and 
B. From Eq. i26\ we see that 



PA(m,i)/Mi 

{A ' B ^ - T X ^ 

A J2 PB (m l )/M t 
i=i 



Vk(AB) ■ 



(27) 



Comparing this expression with Eq. (23) in KG we may iden- 
tify the weights that were introduced in an ad hoc way, in 
order to mimic the effects of large-scale mixing and dilution of 
the SN material. Eq. tells us that & = 1/Mj. Although 
it was shown in KG that the chemical abundance patterns in 
diagrams relating the ratios of two heavier elements (denoted 
A I A diagrams) are insensitive to the weights, we now have a 
theory for calculating these weights. 

The total metallicity distribution fz of low-mass stars, with 
Z = (j4/H), may formally be expressed as a sum of integrals 
using the formalism in KG by replacing their Eq. (15) with our 
new expression for the abundance Uk(A) - However, fz can also 
be found directly by Monte-Carlo simulations, a method used 
in this study. 

For each number k = 1, ...,n we generate Wf. x 10 6 
low-mass stars, where each star has a surface abundance 
of A relative to H determined by Eq. J26i . The weight 
Wk — nk/ 52fc=i n k is the relative number of stars enriched by 
k SNe, where nu is given by Eq. d!6i . The value 10 6 is the 
average number of realizations per number of enriching SNe. 
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Table 1. Input parameters used in application of the theory 



Parameter 


Symbol 


Value 


SN fraction a 


osn 


0.002 


T nw-iruiw star frartion a '^ ) 


alms 


0.835 


Typical density (Halo) c 


T> 


2.05 xl0~ 25 gem" 3 


Density normalization 


Po 


7.18 xl0~ 25 gem -3 


Average SN rate (Halo) d 


Ti 


0.25 kpc" 3 Myr" 1 


SN rate normalization 


Uo 


1.0 kpc~ 3 Myr -1 


Expansion rate c 




7xl0~ 4 kpc 2 Myr -1 



a IMF taken from Kroupa et al. (1991) 
b Stellar lifetimes taken from Portinari et al. (1998) 
Corresponds to a particle density of 0.1 cm -3 
d Estimated from Fig. [4] 

"Expansion rate estimated from Bateman & Larson (1993) 



Thus, the total number of generated stars is n x 10 6 , where 
the maximum number n of enriching SNe is chosen to be large 
enough such that 99.9% of all the low-mass stars enriched by 
n SNe have a metallicity [-A/H] > —2.5. The stellar masses 
rrii are randomly picked from the probability density function 
fm = 4>{in), i-e-, the adopted IMF, and the mixing masses Mi 
are picked from the density functions $M k derived in Eq. \22\ . 
Note that for the stars enriched by k SNe, all k mixing masses 
should be picked from the fcth density function (see Fig. 
All stars enriched by k SNe with the same metallicity are now 
binned together to form, after normalization, the partial proba- 
bility density function fk(Z). Finally, the total probability den- 
sity function is given by 

a 

f z (Z) = J2^kX fk(Z), (28) 

k=l 

where the contribution of /& to fz is given by Wk, the relative 
number of stars enriched by fc SNe (cf. Eq. (22) in KG). 

4. Application of the theory 

In the present theory the chemical enrichment of a metal- 
poor system is completely determined by the five parameter 
functions V m i X (t), ip(t), p(t), 4>(m), andp^rn), whereof three 
(V m i x (t), (j>(m),pA(iri)) are more or less universal (although 
not necessarily known), while the SFR (ip(t)) and the density 
evolution (p(t)) are specific for each system. The intention in 
the following application of the theory is not to present a full 
parameter study, but merely to investigate the effects caused 
by changing the SFR and density evolution. First, however, the 
choice of the other parameters will be commented on. 



chemical enrichment 7 

0.8 1 , , , , , , , , , , , , , 1 
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Time (Gyrs) 

Fig. 4. The rate of Type II SNe in the Galactic halo as a function 
of time. The dots are the result of a chemodynamical simulation 
by Samland et al. (1997), with SN rates given in their Fig. 8. 
A small correction was made for the number of SNe exploding 
in the pre-disk while actually belonging to the halo. The curve, 
given by Eq. ( fJll . is the adopted SN rate in Models B and C. 

into a homogeneous medium assuming no radiative cooling of 
the interior, the radius i?sNR of the remnant is described by 
power-law solutions of the form -Rsnr oc t v , where the value 
i] depends on the phase the remnant is in (Taylor 1950; Sedov 
1959; McKee & Ostriker 1977). However, this initial spread of 
the SN material is confined to a small volume and the chance of 
immediately forming stars in the swept-up gas is probably not 
very large. The successive, large-scale mixing of the cooled SN 
material — caused, e.g., by turbulent motions — must therefore 
be included. Bateman & Larson (1993) estimate the dispersal 
of iron, ascribing the motions in the ISM to different random- 
walk processes where the rms distance i? rms increases as t 1 ^ 2 . 
If so, the mixing volume should on average increase as 

V m Ut) = y r Ls = y (<w0 3/2 kpc 3 , (29) 

where c m i x = | u?r< is a combined "cross-sectional ex- 
pansion rate" for the processes responsible for the bulk dis- 
tribution of heavy elements. These are expected to be cloud 
motions and turbulence in the diffuse medium, assuming the 
spreading in the diluted hot ionized medium can be ne- 
glected (see Bateman & Larson 1993). The Vi in the sum de- 
notes the mean velocity in each phase of the gas, and r, is 
the time between scattering events. From the data given in 
Bateman & Larson (1993) we find an overall expansion rate 
of cr mix = 7x 1CT 4 kpc 2 Myr" 1 (see Table[TJ. Neglecting the 
relatively small initial expansion of the SN remnant itself, an 
approximate estimate of the total mixing volume growth is 
given by Eq. (1291 . where t, measured in Myr, denotes the time 
after the SN explosion. 



4. 1. The mixing volume growth, V mix (t) 

The different evolutionary phases of an SN remnant are well 
understood; and in the ideal case of a blast wave expanding 



4.2. The IMF, <j>{m) 

We adopt the IMF suggested by Pagel (1997). It is a composite 
of the IMF determined by Kroupa et al. (1991) for stars below 
m = 1.0 Mq and two power-law segments approximating the 
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Table 2. Integration times 





Model A 


Model B 


Model C 


n a 


240 


135 


275 




326.4 


200.1 


367.4 


t u Myr c 


1121 


924 


1117 



a Maximum number of enriching SNe (see Sect.|5J 
b Required fi for k — n SNe, according to Eq. 1321 
c Effective upper integration limit for k — n SNe 



20 000 
10 000' 
5 000 

j- 2 000 
1 000 
500 



Model A 
Model B & C 




200 



20 



IMF by Scalo (1986) for stars of higher mass. With the normal- 
ization given by Eq. Q we have 



<j)(m) 



0.506m" - 85 , 0.10 < m < 0.50 

0.253TO" 1 - 85 , 0.50 < m < 1.00 

0.253m- 3 40 , 1.00<m<3.16 

0.113m" 2 - 70 , 3.16 < m < 100. 



(30) 



For this IMF, the fraction of stars that are massive enough 
(> 8 .M©) to eventually explode as core-collapse SNe is 
dSN = 0.002 (see Table Moreover, using the table 
(Z = 0.0004) of stellar lifetimes by Portinari et al. (1998) we 
estimate that (a*) = olms = 0.835, averaged over the first 
billion years. 

4.3. The stellar yield, pA{m) 

Many observational studies use iron (Fe) as a reference ele- 
ment. In order to make a direct comparison easier we shall 
therefore use [Fe/H] as the metallicity indicator. We adopt 
the Fe yields for Population III core-collapse SNe calculated 
by Umeda & Nomoto (2002), assuming that all the SNe in 
the mass interval 13 — 30 M.q have an explosion energy of 
E = 10 51 erg. The fact that the Fe yield is uncertain due to the 
unknown amount of fall-back of 56 Ni onto the neutron star will 
not alter the general conclusions drawn in this paper. 

At later epochs, Type la SNe become an additional source 
of Fe, which is not taken into account in the current model. 
However, observations have shown that up to [Fc/H] ~ —1 
the ISM retained a constant [a/Fe], the a-elements being Mg, 
Si, Ca, and possibly Ti, suggesting that Type II SNe are the 
major source of Fe. Kobayashi et al. (1998) point out that a 
metallicity-dependent rate of Type la SNe could explain this 
relatively late decrease in [a/Fe]. In Sect. 15.51 the model is 
extended to incorporate Fe enrichment by Type la SNe. The 
level of such contamination is found to be small for extremely 
metal-poor stars. 



4.4. The models 

Three models are specified by different sets of the parameters 
and p(t). These parameters will be prescribed, i.e., they 
will not be explicitly calculated from expressions such as Eq. 



40 60 80 

k (no. enriching SNe) 



100 



Fig. 5. The number density of low-mass Galactic halo stars 
as a function of the number fc of enriching SNe. The number 
density of stars decreases with increasing fc, for Model A (filled 
circles), as well as for Model B & C (filled triangles). This 
effect is caused by the fact that the mixing volumes are allowed 
to grow with time. 



<j){m), andpFe(m-) are the same 



<|2J- The parameters V m 
for all three models. 

4.4.1. Model A 



Model A corresponds to the extreme inflow model by Larson 
(1972), which is characterized by a constant SFR/SN rate u 
and a constant gas density p (see Table 0. This means that 
the accretion compensates exactly for the loss of gas in star 
formation. 

The past SFR for the Galactic halo is not very well known. 
Assuming that the rate of SNe in the early Galaxy is compa- 
rable to the present rate of 1 SN every 70 years, we obtain a 
rate of u$n = 0.05 kpc -3 Myr -1 , averaged over the spheroid 
of 40 kpc radius. This is probably a conservative estimate. In 
their chemodynamical simulation of the Galaxy, Samland et al. 
(1997) deduced the rate of Type II SNe in the Galactic halo 
as a function of time (see Fig. |4}. Using their data, we infer 
an SN rate of ugN = 0.25 kpc -3 Myr _1 , averaged over the 
first billion years. This is the adopted rate. Moreover, the aver- 
age particle density is taken to be 0.1 cm -3 , which should be 
a reasonable value for the interstellar gas density of the inner 
Galactic halo, where most of the stars presumably were formed. 

4.4.2. Model B 

Model B is an extension of Model A with an SFR that is al- 
lowed to vary with time. Recalling that ip(t) = usN(0/ a SN, 
we shall assume an SN rate of the form 



«snW = mo 



(t/1000) 
1 + (t/1000) 4 



kpc 3 Myr" 1 , 



(31) 



where t is in Myr and the normalization uq is given in Tabled 
The choice of function is arbitrary and is only used to demon- 
strate the effect of a time-dependent SFR. However, as shown 
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Model B 
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Model C 


/ 
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/ 




x 10 [Fe/H] 



-4.5 -4.0 -3.5 -3.0 -2.5 

[Fe/H] 



Fig. 6. Relative number densities fz of stars as a function of 
metallicity, measured by [Fe/H], for different parameters ip and 
p. See Sect. l4.4l for details. The Fe yields are taken from Umeda 
& Nomoto (2002). 

in Fig. 0] this function resembles the evolution of the SN rate 
given by Samland et al. (1997), with a sharp rise at early times 
followed by a slow decline. As in Model A, the density is given 

by p = p. 
4.4.3. Model C 

In order to investigate the effects of a time-dependent den- 
sity we assume that p oc ip. Thus, p(t) has the same form 
as the expression for usn(£) m Eq. (13 II . with a normalization 
po = 7.18 x 10~ 25 g cm ~ 3 , which corresponds to a maximum 
particle density of 0.2 cm~ 3 at t = 760 Myr. As in Model B, 
the SN rate is given by Eq. ( THV 

5. Results and discussion 



5. 1. Integration times 

Since the probability density function w\sm decreases very fast 
shortly after it has reached its maximum, the expressions for rik 
and fM k do not have to be integrated up to t — tq. An effective 
upper integration limit t u may be introduced. This corresponds 
to an upper limit on p such that p u = p(t u ). The expression 
for p u is chosen to be 



fi u = (ji) + — cr = — Vn+T + n + 1, (32) 

where (p) = n + 1 is the mean and <r = \Jn + 1 is the stan- 
dard deviation of wism for k — n, i.e., the maximum number 
of SNe. With a factor 11/2, the integral of iuism up to p u for 
n = 100 differs only one part in a million from the total in- 
tegral (= 1). Table [2] shows the upper integration limits for 
all three models. Evidently, the population of extremely metal- 
poor stars with a metallicity [Fe/H] < —2.5 is formed after 
~ 1 Gyr, irrespective of the choice of model. 



Table 3. Results 



Model A Model B Model C 



rim kpc 3 


10 250 


7 420 


7 420 


am 


2.0xl0~ 2 


1.4xl0~ 2 


l.OxlO" 2 


7V£ oIII (-4) a kpc- 3 


98 


772 


105 


noIII a 

a Z<-4 


2.2xl0~ 4 


1.5xl0~ 3 


1.5xl0" 4 


F\(-2.5) b 


0.987 


0.996 


0.958 


k d at [Fe/H] = -2.5 C 


142 


74 


158 



a Number/fraction of stars below [Fe/H] = —4, excluding Pop. Ill 
b Fraction of stars enriched by 1 SN below [Fe/H] — —2.5 
c fk d is the dominant contributor to fz at [Fe/H] — —2.5 



5.2. The metallicity distribution 

Before discussing the metallicity distributions we comment 
briefly on the /c-dependence of n^, the number density of (still- 
existing) low-mass stars. As shown in Fig.|5J nk decreases with 
increasing number of enriching SNe. This decrease is present in 
all models. In fact, unless the mixing-volume growth depends 
on the gas density of the ISM (er m ix is assumed to be constant 
in the present study), nk will always decrease with increasing 
k. This result can be ascribed to the properties of M, defined 
in Eq. ( fT4l >. which is a decreasing function of p (and therefore 
time). 



Our resulting metallicity distributions for the three models 
are shown in Fig. [6] with [Fe/H] as metallicity indicator. The 
number density of stars drops drastically below [Fe/H]~ —4. 
This is the low-metallicity tail of stars enriched by a single 
SN. At [Fe/H]> —4, the number density of stars in Models 
A and B approach an exponential increase oc lOFc/H]. Due to 
an initially increasing density (i.e., infall), Model C shows a 
suppressed metal-poor end, similar to the effect caused by in- 
fall in homogeneous models. At [Fe/H]= —2.5, the dominant 
contributing partial density function fk d consists of stars en- 
riched by as many as k^ = 74 — 158 SNe (see Table^i. This 
is in accordance with the estimate by Nissen et al. (1994), who 
concluded that metal-poor Galactic halo stars must have been 
enriched by at least 20 SNe due to the relatively small observed 
star-to-star scatter in the abundance ratios. The recent work by 
Carretta et al. (2002) and Cayrel et al. (2004) verifies a small 
scatter in many abundance ratios down to [Fe/H] = —4. As 
these authors point out, this suggests that even the most metal- 
poor stars observed today may be enriched by more than one 
SN. Figs.0and[8]show that the model stars are indeed enriched 
by several SNe, at least down to [Fe/H] = — 3.5 (cf. Argast et 
al. 2000). This result contradicts the conclusion that most stars 
with [Fe/H] < —2.5 have been enriched by individual SNe 
(e.g., Shigeyama & Tsujimoto 1998). In fact, we cannot be sure 
that a star has been enriched by a single SN unless we observe 
it below the cut-off (see e.g., Fig.0, which in these models is 
located at [Fe/H]~ —4. Note however, that the metallicity dis- 
tributions fk for low k are very broad and about one percent of 
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Fig. 7. The number density of extremely metal-poor stars in the 
Galactic halo as predicted by Model A. The full line is the total 
number density fz, while the dashed lines are the partial num- 
ber densities f k for k = 1,2, 10, and 100 enriching SNe (cf. 
Fig-E) Note the long tail of stars enriched by one SN extending 
beyond [Fc/H] = -2.5. 



1.2 



S0.4 



0.0 



I 



a [Fe/H] =-4.0 



012345678 
k (no. enriching SNe) 




1 2 3 4 5 6 7 
k (no. enriching SNe) 



0.12 



0.018 r 




5 10 15 20 25 30 
k (no. enriching SNe) 



0.000 



50 100 150 200 250 300 
k (no. enriching SNe) 



Fig. 8. The fraction of stars enriched by k SNe at metallicity 
[Fe/H] = -4.0, -3.5, -3.0, and -2.5 as predicted by Model 
A. Note the broad distribution at [Fe/H] = —2.5, peaking at 
kd = 142 (see also Tabled 

Now, in order to compute the fraction of metal-free stars, 
we define the total number density of stars Nz with metallicity 

< Zas 



all the stars enriched by a single SN have presumably a metal- 
licity [Fe/H] > -2.5. This fraction is given by 1 - Fi(-2.5), 
see Table |3 

Admittedly, in the most metal -poor regime the averaging of 
a few contributing SNe cannot be enough to explain the small 
star-to-star scatter, a scatter that also seems to be relatively con- 
stant over a wide range in metallicity. However, we shall argue 
(Paper II) that a small scatter in diagrams relating a ratio of two 
heavier elements to the ratio of an element relative to hydrogen 
(called A/H diagrams in KG) may be a natural consequence of 
selection effects favouring contributions from SNe in a specific 
mass range. 

The simple assumption that the ejected material from k 
SNe mixes within a constant mixing mass (e.g., KG; Fields 
et al. 2002), implying that Fe/H oc k, will underestimate the 
number of enriching SNe at a given metallicity. For example, 
an increase in metallicity by 1.5 dex, from [Fe/H] = —4 to 
[Fe/H] = —2.5, would require a number of k — 30 SNe. This 
is a factor of 2 to 5 less than predicted by the models presented 
here. The reason for this is twofold. Both large-scale mixing of 
the ISM and infall of pristine gas continuously dilute the SN 
material. Thus, in order to reach a certain metallicity in a spe- 
cific region in the ISM, a larger number of SNe, on average, 
have to enrich this region. Note, however, that the distribution 
is very broad, as is evident from Fig.|8] 

5.3. The fraction of Population III stars 

Inserting u>ism(0, t) in the expression for n k we can calculate 
the absolute number density of Population III stars, denoted by 
nk=o = rim. As shown in Table [3] n\\\ = 7420 kpc~ 3 for 
both Models B and C. The higher value for Model A reflects a 
higher average SFR. 



N z (Z) = J2F k (Z)n k , 



(33) 



fe=0 



where F k denotes the distribution function of //. and describes 
the fraction of stars enriched by k SNe with a metallicity < Z. 
It is defined by the integral 



F k {Z) 



f k (Z')dZ'. 



(34) 



The fraction of metal-free stars is then formally given by 



am 



nm 
N z (ooY 



(35) 



Since the integration is terminated at [Fc/H] = —2.5 in the cur- 
rent application, Nz(oo) cannot be calculated directly through 
Eq. d3 3I > - However, making use of the observational finding by 
Carney et al. (1996) that 16% of all the stars have a metallicity 
[Fe/H] < -2.5, we have N z {oo) = 7V z (-2.5)/0.16. Thus, 
am = 0.16nm /Nz{— 2.5). Using this expression, the fraction 
of metal-free stars is estimated to (1 — 2) x 10~ 2 , depending on 
the model (see Table|5Ji- This is close to the prediction by Oey 
(2003), who found a value of am ~ (3 — 4) x 10~ 2 , including 
the effects of additional mixing beyond the local superbubbles. 
This result thus confirms the discrepancy of almost two orders 
of magnitude between the observed and the expected number 
of metal-free stars (Oey 2003). 

If stellar atmospheres are continuously exposed to accre- 
tion of enriched material (e.g., Yoshii et al. 1981), this could, 
at least qualitatively, explain the observed discrepancy. A truly 
metal-free star would then be hidden from us by its polluted at- 
mosphere and convection zone. Another possible explanation is 
that low-mass Population III stars were never formed in the first 
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place (Bond 1981). It has been suggested (e.g., Bromm et al. 
1999; Abel et al. 2002) that the lack of effective cooling agents 
(metals) may inhibit low-mass star formation in the primordial 
ISM. If so, the fraction of stars below, say, [Fe/H] = —4 should 
be given by 



noIII 

a Z<-4 



0.16 



JVz(-4) - mii 
JV z (-2.5)-n m ' 



(36) 



again making use of the observational result by Carney et al. 
(1996). These fractions that exclude Population III stars (see 
Table |5} lie in the interval 10~ 4 — 10~ 3 and are indeed close 
to the present observational value of a^ b <_ 4 = 4 x 10 -4 (cf. 
Oey 2003), which reflects the discovery for HE 0107-5240, the 
most iron-poor star known to date (Christlieb et al. 2002). The 
exact value of _ 4 is as yet, however, uncertain due to small- 
number statistics. 

5.4. The low-metallicity cut-off 

There has been discussion in the past as to whether the ob- 
served cut-off in the number of stars at [Fe/H] ~ —4 (e.g., 
Mc William et al. 1995) is real or whether we should expect a 
continuous distribution down to Z = 0. If there is a cut-off, 
where does it come from? Several authors have argued (e.g., 
Ryan et al. 1996; Shigeyama & Tsujimoto 1998) that a typi- 
cal mixing mass of M ~ 5 x 10 4 A4q should be expected, 
based upon the size of SN remnants when they merge with 
the ambient ISM (Cioffi et al. 1988). This results in a limit- 
ing metallicity of [Fe/H] ~ —2.7, which is at least one order 
of magnitude greater than the limit inferred from observations. 
In fact, this should be regarded as an upper limit (for stars en- 
riched by a single SN), since additional mixing, e.g. due to tur- 
bulent motions in the ISM, would lower the metallicity. On the 
other hand, Oey (2003) estimates a lower limit in the range 
—4.6 < [Fe/H] < —3.6 for the largest superbubbles gener- 
ated by SNe (Oey & Clarke 1997); and Audouze & Silk (1995) 
predict a similar threshold at [Fe/H] ~ —4, assuming that the 
first generation of ~ 10 4 SNe enriched the ISM during a short 
period of time of ~ 2 Myr. 

The statistical approach presented here offers an alternative 
view of the problem. If we allow large-scale mixing of the ISM 
and assume reasonable values for the density and SFR of the 
Galactic halo, the distribution of mixing masses for k = 1 is 
relatively peaked and lies between 10 5 — 10 6 A4q (see Fig.|3}. 
As shown in Fig. [7] this results in a decline in the number of 
stars below [Fc/H] = —4. The definitive value of this thresh- 
old is controlled by the star formation process, which we take 
to be stochastic in nature. Stars enriched by a single SN are 
formed during a finite period of time. During this time they 
trace chemical inhomogeneities of specific sizes, correspond- 
ing to a specific distribution of mixing masses. In fact, contrary 
to the claim that the cut-off is determined by the universal prop- 
erty of SN physics (e.g., Ryan et al. 1996), the present approach 
predicts that a system with, e.g., a different star formation his- 
tory or density evolution produces a different metallicity distri- 
bution with a cut-off determined by the density function fn k=1 , 
the probability of having a mixing mass of size M in regions 
enriched by a single SN. 



5.5. Thermonuclear SNe as an additional Fe source 

The current model primarily describes the chemical enrichment 
of the ISM by the spreading of primary elements produced 
and ejected in core-collapse (Type II) SN explosions. Other 
sources, enriching individual stars in elements with origins dif- 
ferent from core-collapse SNe, e.g., intermediate-mass stars, 
presumably do not have a large impact on the heavy element 
content in the oldest, extremely metal-poor stars. Intermediate- 
mass stars, for example, do not synthesize elements heavier 
than magnesium (excluding the neutron-capture elements). An 
interesting question is, however, whether the large amounts of 
56 Fe ejected in thermonuclear (Type la) SN explosions could 
have an effect on the Fe abundance in extremely metal-poor 
stars, presuming that the rate of Type la SNe is high enough 
at these early epochs. In order to include the chemical enrich- 
ment of the Type la SNe, the current model has to be modified 
slightly. 

5.5.1 . The rate of thermonuclear SNe 

In order to estimate the largest possible contribution of Type la 
SNe to the Fe enrichment in extremely metal-poor stars, we 
adopt a short formation time-scale scenario for the Type la 
SNe. Following Greggio & Renzini (1983) the rate usNia of 
Type la SNe can be estimated by assuming that the SN pre- 
cursors are mass-accreting C/O white dwarfs in close binary 
systems. The formation time-scale of a Type la SN is, in this 
scenario, only determined by the initial mass of the compan- 
ion star. The delay between the onset of star formation and the 
onset of the first Type la SNe is therefore merely ~ 40 Myr, 
corresponding to the lifetime of an 8 M. @ star. 

Assuming that the Galactic SFR follows an exponential 
decay, usNia can be properly normalized by requiring that 
^SNia = usnii at the present epoch (Greggio & Renzini 1983), 
where usnii is the rate of Type II SNe. This means that for, 
e.g., a constant SFR, wsNia — 0.8 x usnii when extrapolated 
to t = tq. The total SN rate replacing Eq. il 1> is given by 
usn = usnu + usNia, where msnii = osn'0- In the discussion 
below, Model A will be used as input the model. However, with 
the inclusion of the Type la SNe, the model will be referred to 
as Model A'. 



5.5.2. The fraction of thermonuclear SNe 

A low-mass star formed at time t has, on average, been en- 
riched by a total number of fi(t) SNe out of which /xsNia(i) 
are SNe of Type la, where /xsNia is given by Eq. (0 
with usn = wsNia (as for ^snii)- Knowing the probability 
/sNia(fc') k) that a star, enriched by a total number of k SNe, 
is enriched by k' Type la SNe (i.e., k — k' Type II SNe), en- 
ables us to determine the total contribution from the Type la 
SNe in detail. The probability of finding a region enriched by 
k' Type la SNe and k — k 1 Type II SNe, respectively, at time t 
is given by Eq. (|4} such that 



WSNIa(fc'.*) 



_ p — MSNIa(t) 



A*SNIa(*)*' /k'l 



(37) 
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a Model A' for stars enriched by k = 100 SNe 



and 

w smi (k - k', t) = e-" w "^fi S mi(t) k - k '/{k - k')\ (38) 
The probability /sNia is then determined by the integral 

TG 

/sNIa(fe', k) (X J W SNIa (fc', t)w S NIl(fc - fc', X 


x a*(£)wisM(M)^(£)dt, (39) 
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[Fe/H] 

Fig. 9. The effect of thermonuclear (Type la) SNe on the metal- 
licity distribution of extremely metal-poor stars. The metal- 
poor tail is completely unaffected by Type la SNe. The frac- 
tion of stars contaminated by one or more Type la SNe (dashed 
line) only becomes significant above [Fe/H] ~ —2.8. The frac- 
tion of stars with > 50% of their Fe produced in Type la SNe 
(dash-dotted line) remains small, even at [Fe/H] = —2.5. The 
metal-rich end of Model A' (full line) is depleted in stars, as 
compared to Model A (dotted line), as these stars are pushed 
towards higher metallicities due to the high Fe yield of the 
Type la SNe. 



where a+(t)wi$,yi{k, t)ip(t) may be regarded as a weighting 
function. Given /sNiai the absolute number density of stars 
enriched by k SNe, of which k 1 are SNe of Type la, is 
n k x /sNia(fc', k). For example, out of 274.1 stars per kpc 3 en- 
riched by 100 SNe (i.e., Model A'), 42.5 kpc -3 are enriched 
by 6 Type la SNe (see Table |4}. The majority of the stars 
enriched by 100 SNe are found to be more metal-rich than 
[Fe/H] = —2.5, since the Fe yield of the Type la SNe is on 
the order of 10 times that of the average Fe yield of Type II 
SNe (e.g., Nomoto et al. 1997). Here, we take the yield to be 
p Fc = 7.49 x 10 _1 M Q for all Type la SNe (Iwamoto et al. 
1999, Model W7). 



5.5.3. The mixing mass distribution for thermo- 
nuclear SNe 

Since the onset of Type la SNe is delayed relative to the Type II 
SNe, the mixing masses of the Type la SNe are, on average, 
smaller than the corresponding mixing masses of the Type II 
SNe at a given time t or for a given total number k of enriching 
SNe. Since all mixing volumes are assumed to grow at the same 
speed, the mixing mass distribution for the Type la SNe are, as 
for the Type II SNe, governed by Eq. (I22> or Eq. d25i . however, 
with the SN rate usn(^ — t~v) replaced by the appropriate rate 
usNia- Similarly, the mixing mass distribution for the Type II 
SNe is now determined by the rate usnii- 



5.5.4. The effect of thermonuclear SNe 

Fig. |9] shows that the metallicity distribution is unaffected by 
the presence of Type la SNe below [Fe/H] ~ —3 for this par- 
ticular choice of input model and ugNia- The metal-rich end of 
Model A' is, however, depleted in stars, as compared to Model 
A. These "missing" stars have been enriched by Type la SNe 
to such an extent that they end up with metallicities higher 
than [Fe/H] = —2.5. Although almost all the remaining stars 
at [Fe/H] = —2.5 may be contaminated by debris from Type la 
SNe (about 95% in Model A'), the Fe in the large majority of 
these stars is produced predominantly in the core-collapse SNe 
(see the dash-dotted curve in Fig.|5}. If the rate of Type la SNe 
was significantly higher in this early phase and if the ISM was 
heavily affected by infall of pristine gas (e.g., as in Model C), 
thermonuclear SNe might have been the dominant source of Fe 
even at metallicities as low as [Fe/H] = —2.5. This is, how- 
ever, at variance with the general belief that Type la SNe did 
not become an important source of Fe until [Fe/H] ~ —1, as 
suggested by observations of stars in the Galactic disk. In or- 
der to achieve such a "late" onset of the Type la SNe, a much 
higher SFR, a much lower density, or a much lower mixing ve- 
locity er mix than adopted would be required in the early Galaxy. 
Note, however, that the low-metallicity cut-off would then be 
displaced accordingly. Alternatively, no (or very few) Type la 
SN explosions occurred in the early Galaxy. In fact, there are 
studies indicating that the formation rate of Type la SN pre- 
cursor systems is metallicity-dependent and that this rate may 
be negligible at metallicities below [Fc/H] = — 1 (Kobayashi 
et al. 1998). In conclusion, extremely metal-poor stars in the 
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[Fe/H] [Mg/H] 

Fig. 10. The difference between high and low psNia/psNii 
yield ratios, a) Metallicity distributions for Fe with (Model A', 
full line) and without (Model A, dashed line) the presence of 
Type la SNe. The onset of the Type la SNe is clearly seen, b) 
The corresponding metallicity distributions, however, with Mg 
as metallicity indicator. The onset of the Type la SNe is virtu- 
ally undetectable due to the low psNia = 8.57 x 10~ 3 A4q of 
Mg (Iwamoto et al. 1999, Model W7). The dotted line denotes 
the metallicity distribution for a 68 times higher psNia of Mg, 
yielding a psNia/psNii ratio similar to that of Fe. 



Galactic halo should, in general, be affected very little by 
Type la SNe. 

Let us take a closer look at the "missing" stars. As men- 
tioned above, these stars are pushed towards higher metallic- 
ities due to the presence of the Type la SNe. In fact, about 
65% of all stars are missing at [Fe/H] = —2.5 in compari- 
son to Model A. Apparently, when the Type la SNe kick in, 
they seem to kick in hard. The effect is caused by the high 
SNIa/SNII yield ratio of Fe (psNia/psNii ~ 9), which creates 
a hole (it actually looks more like a bump when plotted on a 
linear scale, see Fig. 1 10b ) in the metallicity distribution mea- 
sured by Fe. In contrast to holes and bumps created by, e.g., 
variations in the infall rate of pristine gas, this feature will not 
have a distinct counterpart in distributions for elements with a 
lower SNIa/SNII yield ratio, such as Si (psNia/psNii ~ 2), Sc 
(psNia/psNii ~ 1), or even better Mg (psNia/psNii < 0.1, the 
yield ratios are estimated from data published by Iwamoto et al. 
1999). This effect is illustrated in Fig. [TO] It is clear that the on- 
set of the Type la SNe is imprinted differently in metallicity dis- 
tributions of different elements. An Fe hole/bump feature with 
no counterpart in the distribution for, say, Mg, thus suggests an 
additional contribution of Fe from Type la SNe. Such a feature 
may therefore be used to determine whether Galactic halo stars 
with relatively low [a/Fc] ratios have been enriched by Type la 
SNe. Note, however, that for Mg in particular, the metal-rich 
end of the distribution may be affected by intermediate-mass 
stars. 



5.6. Comparison to Oey (2000) 

Although developed independently, our mathematical descrip- 
tion of the evolution of chemical inhomogeneities in the ISM 
resembles that of Oey (2000); nevertheless, the two approaches 
describe two physically different pictures of the chemical en- 
richment in metal-poor systems. We view the evolution of 
chemical inhomogeneities in terms of mixing volumes gen- 
erated by individual SNe. Oey, on the other hand, uses the 
concept of star-forming regions with a fixed distribution of 
sizes V S (N*) of the star-forming regions, where the number 

of SNe within each star-forming region is estimated from 
the mechanical luminosity function, assumed to be a power 
law with index — j3. The metallicity of each star-forming re- 
gion then depends on its predetermined size and is given by 

x my/V s (N*), where m y is the mean yield of metals per 
SN. Conversely, in our approach the size of each mixing vol- 
ume (or mixing mass) grows with time, and the metallicity 
within each volume is determined by the mass-dependent yield 
from a single SN. This enables us to follow the evolution of 
chemical inhomogeneities caused by individual SNe. 

Oey assumes that star-forming regions are uncorrelated and 
randomly distributed in space, meaning that the binomial dis- 
tribution can be used to describe the probability of finding a 
point in the ISM occupied by k (Oey uses j) overlapping re- 
gions. We assume a spatially random distribution as well, but in 
relation to mixing volumes (i.e., uncorrelated star formation). 
Although we use the Poisson distribution instead of the bino- 
mial distribution, the difference is negligible since the average 
number p of enriching SNe is much smaller than the total num- 
ber of SNe in the system. Mathematically, our parameter p cor- 
responds to Oey's parameter nQ, the number of generations of 
star-forming regions times one volume-filling factor per gen- 
eration. However, the two are physically different in that Oey 
assumes that star formation and SN enrichment occurs in local 
regions. Our description can be modified in this direction by 
replacing our mixing volumes with star-forming regions (see 
Paper III). However, we note that our model is not closed; i.e., 
we allow for the gas density to vary with time. 

5. 7. Inhomogeneous infall: The next step 

By treating the infall as a diffuse flow we tend to underesti- 
mate the amplitude of the chemical inhomogeneities. The ini- 
tial and maximum decrease in the local abundance due to an 
infalling cloud of pristine gas is determined by the density con- 
trast p c /p, where p c is the density of the cloud. This ratio does 
not depend on the metallicity of the ISM. On the other hand, 
the abundance increase due to an exploding SN is measured by 
the ratio between the yield p of the element and the mass of 
the element already present in the ISM within the volume af- 
fected by the SN. This ratio can be written as p/eM^ x , where 
e is the mass fraction of the element and is the total mass 
within the affected volume. Comparing these two ratios we see 
that for mixing masses M^ K < (p/e)/(p c / p), the inhomo- 
geneities in the ISM are dominated by those from the SNe. 
Using this relation with, e.g., pp c = 0.1 Mq, £fo = 10 -6 
(corresponding to an abundance [Fe/H] = —3), and p c /p = 1, 



14 



T. Karlsson: Stochastic chemical enrichment 



we obtain M^ x < 10 5 M@. Thus, the large amplitude inho- 
mogeneities in [Fc/H] due to SNe are relatively small in size 
and, as shown in Fig. [3] the SN material is typically mixed 
with a mass which is an order of magnitude larger before be- 
ing locked up in subsequently formed stars. Conversely, if we 
assume that an infalling cloud of mass M c is mixed within a 
volume of mass M^ iK , the criterion for SN-dominated inho- 
mogeneities is M c < (p/e)/(Af 1 ^ c /M^ lix ). Taking the ratio of 
the mixing masses to be on the order of unity, M c < 10 5 A4q 
for the values of pp c and ep c given above. Again, this is a rela- 
tively small mass, in this case for the infalling clouds. Although 
the chemical inhomogeneities are dispersed in the mixing pro- 
cess as well, they are probably affected by heterogeneous infall, 
increasing the star-to-star scatter in abundances relative to hy- 
drogen and displacing the cut-off towards lower metallicities. 
Since a spatially varying density pit, x) would have a similar 
effect on the development of chemical inhomogeneities as in- 
falling gas clouds, the treatment presented in this study gives a 
lower limit to the chemical abundance scatter in stars. However, 
relative abundance ratios of two heavier elements are not in 
themselves affected by how much the gas is diluted. Diagrams 
relating abundance ratios of heavier elements, so-called A/ A 
diagrams, are therefore less sensitive to the mixing uncertain- 
ties than A/H diagrams, in which the abundance of one ele- 
ment is measured relative to hydrogen. 

We note that the present theory may be naturally extended 
to allow for the stochastic infall of gas clouds, in an approach 
similar to that used here for the description of chemical inho- 
mogeneities caused by SNe. 

6. Summary 

The Poisson distribution e - ^' fi(t) k /k\ is taken to describe 
the development chemical inhomogeneities in the ISM caused 
by randomly distributed SNe. If the time-like parameter p(t) 
denotes the average number of SNe that have enriched a ran- 
dom volume element in space at time t, the Poisson distribu- 
tion gives the probability of finding a point in space enriched 
by k SNe at time t. Defining the mixing volume V m i x (t) to be 
the volume at time t affected by material ejected in a single 
SN explosion allows us to account for large-scale mixing by, 
e.g., turbulent motions. This mixing volume is then used in the 
calculation of p(t). Infall is treated as a homogeneous, diffuse 
flow onto the system, affecting the overall gas density pit). 

An analytical expression for the number density of stars 
enriched by a given number k of SNe is derived. It is shown 
that the ratio ip/jft, where ip is the SFR and ^ is the rate of 
change of p, is a measure of n&. An expression for the prob- 
ability density function fj^ k of mixing masses is also given. 
This function describes the probability of how much the ejected 
material from a SN is diluted (i.e., mixed with a certain mix- 
ing mass M), before being locked up in a subsequently formed 
star. This distribution depends on the number k of SNe that 
have enriched the star. 

These results are combined to calculate the metallicity dis- 
tribution of low-mass stars in the Galactic halo. Three mod- 
els with different values for the parameters and pit) are 
discussed (see Sect. 14.4k They all show a steep decline in the 



number of stars below [Fe/H] ~ —4, in agreement with ob- 
servations. This cut-off is a natural consequence of the lower 
probability of finding a large mixing volume enriched by a sin- 
gle SN, producing a small Fe yield. The cut-off is thus deter- 
mined by the tail of fM k=1 (together with the yield). Fig.0and 
Fig. [8] indicate that above the cut-off there is a high probabil- 
ity of finding a star enriched by several SNe, as concluded by 
Nissen et al. (1994) and more recently also suggested, e.g., by 
Cayrel et al. (2004). In fact, only the stars below the cut-off 
have a high probability (~ 1) of being enriched by individual 
SNe. Note, however, that about one percent of all the stars en- 
riched by one single SN have a metallicity [Fc/H] > — 2.5 (see 
Table|3}. 

The absolute number density nni of metal-free 
Population III stars is directly given by Eq. 0, or alter- 
natively Eq. for k = 0. The fraction of metal-free stars 
is estimated to be (1 — 2) x 10~ 2 , which is approximately 
two orders of magnitude higher than the current observational 
value of 4 x 10~ 4 (cf. Oey 2003). This discrepancy was also 
found by Oey (2003). However, if the formation of low-mass 
Population III stars was inhibited due to inefficient cooling of 
the primordial ISM, interrupting further fragmentation (e.g., 
Bromm et al. 1999), the fraction of stars below [Fe/H] = —4 
is brought down to 10 -4 — 10~ 3 , in agreement with the current 
observational estimate. 

Even in the short formation time-scale scenario, and for 
a relatively high rate usisria, the contaminating effect of ther- 
monuclear SNe is found to be insignificant over the entire ex- 
tremely metal-poor regime with a possible exception at the very 
metal-rich end. However, although the majority of stars may be 
contaminated by Type la SNe at [Fe/H] ~ —2.5, the fraction 
of stars with more than 50% of their Fe produced in the Type la 
SNe remains small. Stars that are heavily polluted by Type la 
SNe are pushed towards higher metallicities, thereby creating 
a hole/bump in the metallicity distribution. This effect is pro- 
nounced due to the high yield ratio psNia/fSNii of Fe. Such a 
feature would be virtually absent in distributions for elements 
with much lower psNia/psNii, e.g., for Mg. Hence, the exis- 
tence of an Fe hole/bump feature with no counterpart in, e.g., 
Mg would indicate the presence of a subpopulation of Galactic 
halo stars that have been enriched by Type la SNe. 

Finally, we note that an inhomogeneous infall in the form 
of gas clouds, pristine as well as pre-enriched, would increase 
the amplitude of the chemical inhomogeneities and most proba- 
bly increase the star-to-star scatter in the abundances relative to 
hydrogen. The effect of infalling clouds can be treated stochas- 
tically and should be possible to incorporate in an extended 
version of the present theory. 
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